Thalamocortical Connections Mature along a Sensorimotor-Association Axis of Cortical Developmental Heterochronicty

Read in Data for Developmental Characterization

Glasser parcel names and mesulam assignments

#Glasser regions with corresponding labels, tract names, and mesulam assignments
glasser.regions <- read.csv("/cbica/projects/thalamocortical_development/Maps/parcellations/surface/glasser360_regionlist.csv") #parcel name, label name 
glasser.regions$tract <- paste0("thalamus_", glasser.regions$orig_parcelname) %>% gsub("_ROI", "_autotrack", .) %>% gsub("-", "_", .) #create tract variable with format thalamus_$hemi_$region_autotrack, no dashes -

#Glasser regions assigned to mesulam zones
glasser.assignments <- read.csv("/cbica/projects/thalamocortical_development/Maps/parcellations/surface/glasser360-mesulam_economo_yeo-assignments.csv")
glasser.assignments <- merge(glasser.assignments, glasser.regions, by = "label", sort = F)
glasser.assignments <- glasser.assignments %>% select(label, tract, orig_parcelname, mesulam_assignment)

S-A axis

S.A.axis.cifti <- read_cifti("/cbica/projects/thalamocortical_development//Maps/S-A_ArchetypalAxis/FSLRVertex/SensorimotorAssociation_Axis_parcellated/SensorimotorAssociation.Axis.Glasser360.pscalar.nii") #S-A_ArchetypalAxis repo
S.A.axis <- data.frame(SA.axis = rank(S.A.axis.cifti$data), orig_parcelname = names(S.A.axis.cifti$Parcel))
S.A.axis <- merge(S.A.axis, glasser.assignments, by = c("orig_parcelname"), sort = F)
S.A.axis$SA.axis.bin <- as.numeric(cut2(S.A.axis$SA.axis, g = 5)) #divide the S-A axis into 5 ranked bins, 72 regions per bin

S.A.axis %>% group_by(SA.axis.bin) %>% do(num.regions = length(.$SA.axis), mean.rank = mean(.$SA.axis)) %>% unnest(cols=c("num.regions", "mean.rank"))
## # A tibble: 5 × 3
##   SA.axis.bin num.regions mean.rank
##         <dbl>       <int>     <dbl>
## 1           1          72      36.5
## 2           2          72     108. 
## 3           3          72     180. 
## 4           4          72     252. 
## 5           5          72     324.

NeuroSynth cognitive atlas term maps

#Neurosynth z-score maps generated by nimare. Nimare uses multilevel kernel density analysis- Chi-square analysis + FDR-correction to use the same procedure as Neurosynth
neurosynth.terms <- read.csv("/cbica/projects/thalamocortical_development/Maps/neurosynth/atl-glasser360_neurosynth.csv") %>% select(-regionName) %>% select(-timing) #neurosynth term meta-analytic scores for 123 terms present in both NeuroSynth and the Cognitive Atlas, ordered in lh --> rh glasser order
neurosynth.termlist <- names(neurosynth.terms) #list of terms (cue willy wonka line "dont you want to know our names?" "cant imagine why it wouldnt matter")
neurosynth.terms$label[1:180] <- glasser.regions$label[181:360] #lh first
neurosynth.terms$label[181:360] <- glasser.regions$label[1:180] #rh

Spatial permutation (spin) test parcel rotation matrix

perm.id.full <- readRDS("/cbica/projects/thalamocortical_development/software/rotate_parcellation/glasser_sphericalrotations_N10000.rds") #10,000 spatial null spins
spin.parcels <- glasser.regions #order of complete set of glasser parcels for spinning

Dataset-specific diffusion measure spreadsheets

FA.glasser.pnc <- read.csv("/cbica/projects/thalamocortical_development/thalamocortical_structuralconnectivity/individual/PNC/PNC_FA_finalsample_primary.csv") #generated by sample_construction/PNC/tractmeasures_dfs_PNC.R

FA.glasser.hcpd <- read.csv("/cbica/projects/thalamocortical_development/thalamocortical_structuralconnectivity/individual/HCPD/HCPD_FA_finalsample_primary.csv") #generated by /sample_construction/HCPD/tractmeasures_dfs_HCPD.R

Dataset-specific tract lists

tractlist.pnc <- read.csv("/cbica/projects/thalamocortical_development/thalamocortical_structuralconnectivity/individual/PNC/PNC_tractlist_primary.csv") #generated by tract_measures/tractlists/thalamocortical_tractlists.R

tractlist.hcpd <- read.csv("/cbica/projects/thalamocortical_development/thalamocortical_structuralconnectivity/individual/HCPD/HCPD_tractlist_primary.csv") #generated by tract_measures/tractlists/thalamocortical_tractlists.R

Developmental statistics

setwd("/cbica/projects/thalamocortical_development/thalamocortical_results/development_results/") #gam output dir

#list all FA development files in development results directory for accessing. Files were generated by gam_functions/fit_ageGams.R
files <- list.files(getwd(), pattern='FA', ignore.case=T, full.names = F) 
#files <- files[!str_detect(files,pattern="temporalSAcorr")] #except these, which we will read in later because they are very large

filenames <- c()
#generate variable names to assign files data to 
for(name in files){
  Rname <- gsub('^.{12}|.{4}$', '', name) #remove ".csv" from end of filename and "development_" from start of filename
  filenames <- append(filenames, Rname) #save filenames into a character vector 
}

#read in files and assign to variables
for(i in 1:length(filenames)){
  
  Rfilename <- sprintf("%s",filenames[i]) #save filename as an individual string
  
  if(grepl("gameffects", Rfilename)) {
    x <- read.csv(files[i], header = TRUE) 
    x <- merge(x, glasser.assignments, by = "tract")
    assign(Rfilename, x) 
  }
  
  if(!grepl("temporal", Rfilename) & !grepl("gameffects", Rfilename)) {
    x <- read.csv(files[i], header = TRUE) 
    x <- merge(x, S.A.axis, by = "tract")
    assign(Rfilename, x) 
  }
  
  if(grepl("temporal", Rfilename)){
    x <- read.csv(files[i], header = FALSE) 
    assign(Rfilename, x) 
  }
}
rm(x)

A Spectrum of Thalamocortical Developmental Trajectories

Cortex-wide thalamic connectivity developmental profiles

Number of significant developmental effects

#Function to calculate the number and percent of thalamocortical connections showing a significant developmental change in a given measure
sig.effects <- function(measure, atlas, dataset){
  ageeffects.df <- get(sprintf("gameffects_%s_%s_%s", measure, atlas, dataset)) 
  ageeffects.df$significant <- p.adjust(ageeffects.df$Anova.smooth.pvalue, method = c("fdr")) < 0.05 #fdr-corrected significant connections
  sigeffect.totaln <- sum(ageeffects.df$significant) #total number of significant connections 
  sigeffect.percent <- round(sigeffect.totaln/length(ageeffects.df$significant)*100, 2) #percent of significant connections
  print(sprintf("In %s, %s thalamocortical connections (%s percent) show significant age-related changes in %s", dataset, sigeffect.totaln, sigeffect.percent, measure))
}

PNC

sig.effects(measure = "FA", atlas = "glasser", dataset = "pnc")
## [1] "In pnc, 201 thalamocortical connections (90.13 percent) show significant age-related changes in FA"

HCPD

sig.effects(measure = "FA", atlas = "glasser", dataset = "hcpd")
## [1] "In hcpd, 174 thalamocortical connections (75.65 percent) show significant age-related changes in FA"

Thalamocortical connection GAM smooth functions

PNC

smoothcentered_FA_glasser_pnc.plot <- merge((smoothcentered_FA_glasser_pnc %>% filter(grepl("thalamus_R", tract))), (gameffects_FA_glasser_pnc %>% select(tract, GAM.smooth.partialR2)), by = "tract") #zero-centered developmental smooth functions for RH connections + tract-specific partial R2 for plotting

ggplot(smoothcentered_FA_glasser_pnc.plot, aes(x = age, y = est, group = tract, color = GAM.smooth.partialR2)) + 
  geom_line(linewidth = .3, alpha = .8) + 
  scale_color_gradientn(colors = c("#FEC22F", "#F59A72", "#9c3a80", "#672975","#323280"), limits = c(0, 0.1), oob = squish) +
  theme_minimal() +
  labs(x="\nAge", y="Connection FA\n") +
  scale_x_continuous(breaks=c(8, 10, 12, 14, 16, 18, 20, 22), expand = c(0,.45)) +
  theme(
  legend.position = "none", 
  panel.grid.major = element_blank(),
  panel.grid.minor = element_blank(),
  axis.text = element_text(size = 5, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2))

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure4/Developmental_Trajectories_PNC.pdf", device = "pdf", dpi = 500, width = 2.02, height = 1.65)

HCPD

smoothcentered_FA_glasser_hcpd.plot <- merge((smoothcentered_FA_glasser_hcpd %>% filter(grepl("thalamus_R", tract))), (gameffects_FA_glasser_hcpd %>% select(tract, GAM.smooth.partialR2)), by = "tract")

ggplot(smoothcentered_FA_glasser_hcpd.plot, aes(x = age, y = est, group = tract, color = GAM.smooth.partialR2)) + 
  geom_line(linewidth = .3, alpha = .8) + 
  scale_color_gradientn(colors = c("#FEC22F", "#F59A72", "#9c3a80", "#672975","#323280"), limits = c(-0.02, 0.075), oob = squish) +
  theme_minimal() +
  labs(x="\nAge", y="Connection FA\n") +
  scale_x_continuous(breaks=c(8, 10, 12, 14, 16, 18, 20, 22), expand = c(0,.45)) +
  scale_y_continuous(limits = c(-0.02, 0.0138)) +
  theme(
  legend.position = "none", 
  panel.grid.major = element_blank(),
  panel.grid.minor = element_blank(),
  axis.text = element_text(size = 5, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2))

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure4/Developmental_Trajectories_HCPD.pdf", device = "pdf", dpi = 500, width = 2.02 , height = 1.65)

Region-specific GAM smooths

#Function for connection-specific GAM spline + derivative plots
connection_trajectory_plot <- function(measure, atlas, dataset, tract_name, display_color, y_ticks){
  
  #Read in participant-level data and model-level fitted values and derivatives 
  participant.data.df <- get(sprintf("%s.%s.%s", measure, atlas, dataset))
  fittedvalues.df <- get(sprintf("fittedvalues_%s_%s_%s", measure, atlas, dataset)) %>% filter(., tract == tract_name)
  derivatives.df <- get(sprintf("derivatives_%s_%s_%s", measure, atlas, dataset)) %>% filter(., tract == tract_name)
  derivatives.df$significant.derivative[derivatives.df$significant.derivative == 0] <- NA

  #Connection spline plot with participant-level data
  trajectory.plot <- ggplot(data = participant.data.df, aes(x = age, y = get(tract_name))) +
    geom_point(color = alpha(display_color, 0.5), size = .1) +
    geom_line(data = fittedvalues.df, aes(x = age, y = fitted), color = display_color, linewidth = 1) +
    geom_ribbon(data = fittedvalues.df, aes(x = age, y = fitted,  ymin = lower, ymax = upper), alpha = .8, linetype = 0, fill = display_color) +
    labs(x="\nAge", y=sprintf("%s %s: %s\n", tract_name, measure, toupper(dataset))) +
    theme_classic() +
    scale_x_continuous(breaks = c(8, 12, 16, 20), limits = c(8,23), expand = c(0.025,.45)) +
    scale_y_continuous(breaks = y_ticks) + 
    theme(
        legend.position = "none", 
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.text = element_text(size = 5, family = "Arial", color = c("black")),
        axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
        axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
        axis.line = element_line(linewidth = .2), 
        axis.ticks = element_line(linewidth = .2))
  
  #Significant derivatives (rate of age-related change) plot 
  derivatives.plot <- ggplot(data = derivatives.df) + 
    geom_tile(aes(x = age, y = .1, fill = significant.derivative)) +
    scale_fill_gradient(low = alpha(display_color, 0.2), high = display_color,  na.value = "white") +
    labs(x="\nAge", fill = "Derivative") + 
    theme_classic() +
    scale_x_continuous(breaks=c(8, 12, 16, 20), limits = c(8,23), expand = c(0.025,.45)) +
    theme(
        legend.position = "none",
        axis.title.y = element_blank(),
          axis.text = element_blank(),
          axis.line = element_blank(),
          axis.ticks = element_blank(),
          axis.title = element_blank(),
          legend.text =  element_text(size = 6, family = "Arial", color = c("black"))) 
  
  allplots <- list(trajectory.plot, derivatives.plot)
  tract.plot <- cowplot::plot_grid(rel_heights = c(16, 3.2), plotlist = allplots, align = "v", axis = "lr", ncol = 1)

  return(tract.plot) 
}

Primary motor cortex (4)

Left 4: S-A axis rank = 16, mesulam assignment = idiotypic

connection_trajectory_plot(measure = "FA", atlas = "glasser", dataset = "pnc", tract_name = "thalamus_L_4_autotrack", display_color = "#FEC22F", y_ticks = c(0.4, 0.45, 0.5, 0.55))

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure4/PrimaryMotor_PNC.pdf", device = "pdf", dpi = 500, width = 1.95, height = 1.8)
gameffects_FA_glasser_pnc %>% filter(tract == "thalamus_L_4_autotrack") %>% select(smooth.increase.offset)
##   smooth.increase.offset
## 1               16.06784
connection_trajectory_plot(measure = "FA", atlas = "glasser", dataset = "hcpd", tract_name = "thalamus_L_4_autotrack", display_color = "#FEC22F", y_ticks = c(0.35, 0.4, 0.45, 0.5))

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure4/PrimaryMotor_HCPD.pdf", device = "pdf", dpi = 500, width = 1.95, height = 1.8)
gameffects_FA_glasser_hcpd %>% filter(tract == "thalamus_L_4_autotrack") %>% select(smooth.increase.offset)
##   smooth.increase.offset
## 1               15.38233

Lateral parietal (PF)

Left PF: S-A axis rank = 248, mesulam assignment = heteromodal

connection_trajectory_plot(measure = "FA", atlas = "glasser", dataset = "pnc", tract_name = "thalamus_L_PF_autotrack", display_color = "#672975", y_ticks = c(0.3, 0.35, 0.4, 0.45))

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure4/LateralParietal_PNC.pdf", device = "pdf", dpi = 500, width = 1.95, height = 1.8)
gameffects_FA_glasser_pnc %>% filter(tract == "thalamus_L_PF_autotrack") %>% select(smooth.increase.offset)
##   smooth.increase.offset
## 1                18.4531
connection_trajectory_plot(measure = "FA", atlas = "glasser", dataset = "hcpd", tract_name = "thalamus_L_PF_autotrack", display_color = "#672975", y_ticks = c(0.3, 0.35, 0.4, 0.45))

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure4/LateralParietal_HCPD.pdf", device = "pdf", dpi = 500, width = 1.95, height = 1.8)
gameffects_FA_glasser_hcpd %>% filter(tract == "thalamus_L_PF_autotrack") %>% select(smooth.increase.offset)
##   smooth.increase.offset
## 1               16.42504

Lateral prefrontal (8BL)

Left 8BL: S-A axis rank = 342, mesulam assignment = heteromodal

connection_trajectory_plot(measure = "FA", atlas = "glasser", dataset = "pnc", tract_name = "thalamus_L_8BL_autotrack", display_color = "#323280", y_ticks = c(0.3, 0.35, 0.4, 0.45))

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure4/LateralFrontal_PNC.pdf", device = "pdf", dpi = 500, width = 1.95, height = 1.8)
gameffects_FA_glasser_pnc %>% filter(tract == "thalamus_L_8BL_autotrack") %>% select(smooth.increase.offset)
##   smooth.increase.offset
## 1               21.21106
connection_trajectory_plot(measure = "FA", atlas = "glasser", dataset = "hcpd", tract_name = "thalamus_L_8BL_autotrack", display_color = "#323280", y_ticks = c(0.3, 0.35, 0.4, 0.45))

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure4/LateralFrontal_HCPD.pdf", device = "pdf", dpi = 500, width = 1.95, height = 1.8)
gameffects_FA_glasser_hcpd %>% filter(tract == "thalamus_L_8BL_autotrack") %>% select(smooth.increase.offset)
##   smooth.increase.offset
## 1               21.91667

Maturational patterns are highly reproducible

Correlation of thalamocortical connection development effect sizes between datasets

Pearson’s R

gameffects.merged <- merge(gameffects_FA_glasser_pnc, gameffects_FA_glasser_hcpd, by = c("tract", "label", "orig_parcelname", "mesulam_assignment"), suffixes = c("_pnc", "_hcpd"))

cor.test(gameffects.merged$GAM.smooth.partialR2_pnc, gameffects.merged$GAM.smooth.partialR2_hcpd)
## 
##  Pearson's product-moment correlation
## 
## data:  gameffects.merged$GAM.smooth.partialR2_pnc and gameffects.merged$GAM.smooth.partialR2_hcpd
## t = 15.505, df = 219, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.6538415 0.7808211
## sample estimates:
##       cor 
## 0.7233925

Spatial permutation (spin) based p-value

spin.data <- left_join(spin.parcels, gameffects.merged, by = c("label", "orig_parcelname", "tract")) #full set of parcel data in rh --> lh order for spinning. spin test null correlations use complete obs only. Each null correlation correspondence statistic is thus computed on a slightly  reduced set of data, akin to a jackknife procedure
perm.sphere.p(x = spin.data$GAM.smooth.partialR2_pnc, y = spin.data$GAM.smooth.partialR2_hcpd, perm.id = perm.id.full, corr.type = "pearson") 
## [1] 0

Correlation plot

ggplot(gameffects.merged, aes(x = GAM.smooth.partialR2_pnc, y = GAM.smooth.partialR2_hcpd)) +
  geom_point(color = c("#FCAB6A"), size = 0.5) +
  geom_smooth(method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) + 
  theme_classic() +
  scale_x_continuous(limits = c(-0.01, 0.17)) +
  scale_y_continuous(limits = c(-0.06, 0.15)) +
  labs(x="\nPNC", y="HCPD\n") +
  theme(
  axis.text = element_text(size = 5, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2)) 

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure4/FA_Ageeffect_correlation_PNCHCPD.pdf", device = "pdf", dpi = 500, width = 1.99, height = 1.58)

Correlation of thalamocortical connection age of maturation between datasets

Pearson’s R

cor.test(gameffects.merged$smooth.increase.offset_pnc, gameffects.merged$smooth.increase.offset_hcpd)
## 
##  Pearson's product-moment correlation
## 
## data:  gameffects.merged$smooth.increase.offset_pnc and gameffects.merged$smooth.increase.offset_hcpd
## t = 7.8687, df = 150, p-value = 6.581e-13
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.4172322 0.6442871
## sample estimates:
##      cor 
## 0.540529

Spatial permutation (spin) based p-value

perm.sphere.p(x = spin.data$smooth.increase.offset_pnc, y = spin.data$smooth.increase.offset_hcpd, perm.id = perm.id.full, corr.type = "pearson") 
## [1] 0.00025

Correlation plot

ggplot(gameffects.merged, aes(x = smooth.increase.offset_pnc, y = smooth.increase.offset_hcpd)) +
  geom_point(color = c("#9c3a80"), size = 0.4) +
  geom_smooth(method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) + 
  theme_classic() +
  scale_x_continuous(limits = c(14, 23), breaks = c(14, 16, 18, 20, 22), expand = c(0.05, 0)) +
  scale_y_continuous(limits = c(14, 23), breaks = c(14, 16, 18, 20, 22), expand = c(0, 1)) +
  labs(x="\nPNC", y="HCPD\n") +
  theme(
  axis.text = element_text(size = 5, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2)) 

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure4/FA_Agematuration_correlation_PNCHCPD.pdf", device = "pdf", dpi = 500, width = 1.94, height = 1.58)

Correlation of developmental magnitude and age of maturation within dataset

PNC

Spearman’s rho

cor.test(gameffects_FA_glasser_pnc$GAM.smooth.partialR2, gameffects_FA_glasser_pnc$smooth.increase.offset, method = c("spearman"), exact = F)
## 
##  Spearman's rank correlation rho
## 
## data:  gameffects_FA_glasser_pnc$GAM.smooth.partialR2 and gameffects_FA_glasser_pnc$smooth.increase.offset
## S = 229591, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.8278022

Spatial permutation (spin) based p-value

spin.data <- left_join(spin.parcels, gameffects_FA_glasser_pnc, by = c("label", "orig_parcelname", "tract")) 
perm.sphere.p(x = spin.data$GAM.smooth.partialR2, y = spin.data$smooth.increase.offset, perm.id = perm.id.full, corr.type = "spearman") 
## [1] 0

Correlation plot

ggplot(gameffects_FA_glasser_pnc, aes(x = GAM.smooth.partialR2, y = smooth.increase.offset)) +
  geom_point(color = c("#FCAB6A"), size = 0.5) +
  geom_smooth(method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) + 
  theme_classic() +
  scale_y_continuous(limits = c(12, 23), breaks = c(12, 14, 16, 18, 20, 22)) +
  labs(x="\nAge effect (partial R2)", y="Age of maturation (years)\n") +
  theme(
  axis.text = element_text(size = 5, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2)) 

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure4/FA_Ageeffect_Ageofmaturation_PNC.pdf", device = "pdf", dpi = 500, width = 1.93, height = 1.58)

HCPD

cor.test(gameffects_FA_glasser_hcpd$GAM.smooth.partialR2, gameffects_FA_glasser_hcpd$smooth.increase.offset, method = c("spearman"), exact = F)
## 
##  Spearman's rank correlation rho
## 
## data:  gameffects_FA_glasser_hcpd$GAM.smooth.partialR2 and gameffects_FA_glasser_hcpd$smooth.increase.offset
## S = 122759, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.8267486

Spatial permutation (spin) based p-value

spin.data <- left_join(spin.parcels, gameffects_FA_glasser_hcpd, by = c("label", "orig_parcelname", "tract")) 
perm.sphere.p(x = spin.data$GAM.smooth.partialR2, y = spin.data$smooth.increase.offset, perm.id = perm.id.full, corr.type = "spearman") 
## [1] 0

Early and Late Maturing Thalamocortical Connections

Early and late maturing areas

#Identify parcels that fall within first and fourth quartiles of the age of maturation variable
matstage.df <- gameffects_FA_glasser_pnc %>% select(label, smooth.increase.offset)
matstage.df$maturation_bracket <- "mid" #set default bracket to mid
matstage.df$maturation_bracket[matstage.df$smooth.increase.offset < quantile(matstage.df$smooth.increase.offset, 0.25, na.rm = T)] <- "early" #set first quartile bracket to early
matstage.df$maturation_bracket[matstage.df$smooth.increase.offset > quantile(matstage.df$smooth.increase.offset, 0.75, na.rm = T)] <- "late" #set fourth quartile bracket to late

matstage.df <- rbind(matstage.df, data.frame(label = "rh_???", smooth.increase.offset = NA, maturation_bracket = "medialwall")) #create a label for medial wall to color it grey
matstage.df <- rbind(matstage.df, data.frame(label = "lh_???", smooth.increase.offset = NA, maturation_bracket = "medialwall")) #create a label for medial wall to color it grey

ggseg(.data = matstage.df, atlas = "glasser", mapping = aes(fill = maturation_bracket, colour=I("black"), size=I(.06))) + scale_fill_manual(values = c("#FEC22F", alpha("#323280", 0.97), "white", "white" ), na.value = "gray90") + theme(legend.position = "none")
## merging atlas and data by 'label'
## Warning: Some data not merged properly. Check for naming errors in data:
##   atlas type  hemi  side  region label smooth.increase.offset maturation_bracket
##   <chr> <chr> <chr> <chr> <chr>  <chr>                  <dbl> <chr>             
## 1 <NA>  <NA>  <NA>  <NA>  <NA>   lh_L…                   16.9 mid

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure4/Maturationalbrackets_PNC.pdf", device = "pdf", dpi = 500, width = 4.4, height = 4.4)

NeuroSynth decoding of age of maturation maps

Compute cognitive term-development map spatial correlations

#Function to calculate the correlation between a neurosynth term z-score map and a thalamocortical development map
developmentmap.neurosynth.decoding <- function(df, dev.metric, term){
  df.allregions <- left_join(spin.parcels, df, by = c("label", "tract", "orig_parcelname")) #full set of parcel data in rh --> lh order 
  term.cor <- cor(df.allregions[,term], df.allregions[,dev.metric], use = "complete.obs", method = c("spearman")) #correlation between neurosynth term map and development metric
  term.results <- data.frame("term" = term, "term.correlation" = term.cor)
  return(term.results)
}

Developmental timing decoding by Neurosynth cognitive terms

PNC

gameffects_neurosynth_pnc <- merge(gameffects_FA_glasser_pnc, neurosynth.terms, by = "label")

devpattern.neurosynth.pnc <- ldply(neurosynth.termlist, function(t){
  developmentmap.neurosynth.decoding(df = gameffects_neurosynth_pnc, dev.metric = "smooth.increase.offset", term = t)}) %>% #neurosynth-neurodevelopment correlations for all terms
  arrange(desc(term.correlation))
devpattern.neurosynth.pnc <- rbind(slice_max(devpattern.neurosynth.pnc, order_by = term.correlation, n = 10), slice_min(devpattern.neurosynth.pnc, order_by = term.correlation, n = 10)) #10 most positively correlated and negatively correlated neurosynth terms with the age of maturation map

devpattern.neurosynth.pnc
##                   term term.correlation
## 1           expectancy        0.4111486
## 2           monitoring        0.3976832
## 3    cognitive_control        0.3815389
## 4            reasoning        0.3296662
## 5             strategy        0.3109975
## 6  response_inhibition        0.3094288
## 7     memory_retrieval        0.2866834
## 8              thought        0.2835596
## 9               memory        0.2775858
## 10              recall        0.2677131
## 11  object_recognition       -0.2935350
## 12   visual_perception       -0.2879558
## 13          morphology       -0.2742436
## 14       consolidation       -0.2560287
## 15                gaze       -0.2521194
## 16        coordination       -0.2340298
## 17           expertise       -0.2327537
## 18          perception       -0.2242119
## 19   facial_expression       -0.2200723
## 20           induction       -0.2038560

HCPD

gameffects_neurosynth_hcpd <- merge(gameffects_FA_glasser_hcpd, neurosynth.terms, by = "label")

devpattern.neurosynth.hcpd <- ldply(neurosynth.termlist, function(t){
  developmentmap.neurosynth.decoding(df = gameffects_neurosynth_hcpd, dev.metric = "smooth.increase.offset", term = t)}) %>%
  arrange(desc(term.correlation))
devpattern.neurosynth.hcpd <- rbind(slice_max(devpattern.neurosynth.hcpd, order_by = term.correlation, n = 10), slice_min(devpattern.neurosynth.hcpd, order_by = term.correlation, n = 10))

devpattern.neurosynth.hcpd
##                      term term.correlation
## 1              expectancy        0.3940717
## 2  sentence_comprehension        0.3852283
## 3                 thought        0.3694649
## 4       cognitive_control        0.3561159
## 5     response_inhibition        0.3480972
## 6              monitoring        0.3193377
## 7      emotion_regulation        0.3186345
## 8                  recall        0.2951983
## 9                language        0.2831427
## 10             inhibition        0.2653346
## 11     object_recognition       -0.3456587
## 12           multisensory       -0.3419028
## 13      visual_perception       -0.3341608
## 14             adaptation       -0.2941774
## 15         mental_imagery       -0.2867927
## 16              induction       -0.2677590
## 17           localization       -0.2631644
## 18               fixation       -0.2611707
## 19                   gaze       -0.2601696
## 20                imagery       -0.2521325

Overlapping terms across datasets

merge(devpattern.neurosynth.pnc, devpattern.neurosynth.hcpd, by = "term") %>% arrange(desc(term.correlation.x)) %>% select(term)
##                   term
## 1           expectancy
## 2           monitoring
## 3    cognitive_control
## 4  response_inhibition
## 5              thought
## 6               recall
## 7            induction
## 8                 gaze
## 9    visual_perception
## 10  object_recognition
ggplot(devpattern.neurosynth.pnc, aes(x = term.correlation, y = term, color = term.correlation)) +
  geom_segment(aes(y = reorder(term, term.correlation), yend = term, x = 0, xend =  
                     term.correlation), color = "#989898", linewidth = 0.2) +
  geom_point(size = 2.2, alpha = 0.85) +
  theme_light() +
  scale_color_gradientn(colors = c("#FEC22F", "#F59A72", "#9c3a80", "#672975","#323280"), limits = c(-0.29, 0.35), oob=squish) +
  labs(x = "\nCognitive Term-Thalamocortical Development Correlation", y = "NeuroSynth Term\n") +
  scale_x_continuous(limits = c(-0.3,0.42), breaks = c(-0.3, -0.2, -0.1, 0, 0.1, 0.2, 0.3, 0.4)) +
  scale_y_discrete(
  labels = c("expectancy"=expression(bold(expectancy)),
  "monitoring"=expression(bold(monitoring)),
  "cognitive_control"=expression(bold(cognitive_control)),
  "response_inhibition"=expression(bold(response_inhibition)),
  "thought"=expression(bold(thought)),
  "recall"=expression(bold(recall)),
  "induction"=expression(bold(induction)),
  "object_recognition"=expression(bold(object_recognition)),
  "visual_perception"=expression(bold(visual_perception)),
  "gaze"=expression(bold(gaze)),
  "induction"=expression(bold(induction), parse=TRUE))) +
   theme(
   legend.position = "none", 
   panel.grid.major.y = element_blank(),
   panel.grid.major.x = element_line(color = c("gray90")),
   axis.ticks.y = element_blank(),
   panel.border = element_blank(),
   axis.text = element_text(size = 5, family = "Arial", color = c("black")),
   axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
   axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
   axis.ticks = element_line(linewidth = .2))

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure4/Neurosynth_Maturationdecoding.pdf", device = "pdf", dpi = 500, width = 3.35, height = 3.07)

Probability of obtaining 10/20 overlapping words given the set of 123 terms

#Function to compute the number of overlapping terms when drawing term sets of size 20 at random for PNC and HCPD nulls
compute_term_overlap <- function(){
term.subset.pnc <- sample(x = neurosynth.termlist, size = 20, replace = FALSE)
term.subset.hcpd <- sample(x = neurosynth.termlist, size = 20, replace = FALSE)
num.overlap <- length(intersect(term.subset.pnc, term.subset.hcpd))
return(num.overlap)}

#Null distribution of overlapping term counts based on 10,000 random draws
set.seed(1)
overlap.number <- replicate(10000, compute_term_overlap()) %>% as.data.frame %>% set_names("overlaps")

sprintf("In %s out of 10,000 runs, the term overlap number was equal to or greater than 10", sum(overlap.number$overlaps >= 10))
## [1] "In 1 out of 10,000 runs, the term overlap number was equal to or greater than 10"
sprintf("The permutation-based p value for these overlap counts is %s", sum(overlap.number$overlaps >= 10) / 10000)
## [1] "The permutation-based p value for these overlap counts is 1e-04"
summary(overlap.number$overlaps)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.000   2.000   3.000   3.257   4.000  10.000
ggplot(overlap.number, aes(x = overlaps)) + 
  geom_histogram(bins = 12, fill = c( "#9c3a80"), color = c("white")) +
  theme_classic() +
  scale_x_continuous(limits = c(-1, 10), breaks = c(0, 2, 4, 6, 8, 10)) +
  labs(x = "\nNeurosynth term overlap", y = "Count\n") +
  geom_vline(xintercept = 10, linewidth = .1) +
  theme(
  legend.position = "none", 
  axis.text = element_text(size = 5, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .1), 
  axis.ticks = element_line(linewidth = .1))

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure4/Neurosynth_overlap_null.pdf", device = "pdf", dpi = 500, width = 1.3, height = 1.1)

Thalamocortical Connectivity Maturation is Organized by the Sensorimotor-Association Axis

gameffects_FA_glasser_pnc <- merge(gameffects_FA_glasser_pnc, (S.A.axis %>% select(tract, SA.axis)), by = "tract")
gameffects_FA_glasser_hcpd <- merge(gameffects_FA_glasser_hcpd, (S.A.axis %>% select(tract, SA.axis)), by = "tract")

Age-resolved developmental change plots

PNC

derivatives_FA_glasser_pnc$significant.derivative[derivatives_FA_glasser_pnc$significant.derivative <= 0] <- NA #replace non-significantly increasing derivatives with NA. (Note here 0 was == not significant)

ggplot() +
  geom_line(data = derivatives_FA_glasser_pnc, aes(x = age, y = SA.axis, group = tract, alpha = significant.derivative, color = SA.axis, linewidth = significant.derivative)) +
  scale_alpha_continuous(range = c(0, 1), breaks = c(0, 0.001, 0.002, 0.003, 0.004, 0.005, 0.006), limits = c(0.0005, 0.006)) +
  scale_linewidth_continuous(range = c(0, 4.5), breaks = c(0, 0.001, 0.002, 0.003, 0.004, 0.005, 0.006),  limits = c(0.0005, 0.006)) +
  scale_color_gradient2(low = "#ffc933", mid = "#e9dcf2", high = "#6f1282", guide = "colorbar", midpoint = 180)  +
  scale_x_continuous(breaks = c(8, 10, 12, 14, 16, 18, 20, 22), expand = c(0, .25)) +
  scale_y_continuous(expand = c(0.03, 0)) +
  theme_classic() +
  ylab("Sensorimotor-Association Axis Rank\n") +
  xlab("\nAge") +
  theme(
   legend.position = "none",
   axis.line = element_line(size = .2),
   axis.text = element_text(size = 5, family = "Arial", color = c("black")),
   axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
   axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
   axis.ticks = element_line(linewidth = .2))

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure5/Derivativeplot_pnc.pdf", device = "pdf", dpi = 500, width = 2.3, height = 3)

HCPD

derivatives_FA_glasser_hcpd$significant.derivative[derivatives_FA_glasser_hcpd$significant.derivative <= 0] <- NA #replace non-significantly increasing derivatives with NA. (Note here 0 was == not significant)

ggplot() +
  geom_line(data = derivatives_FA_glasser_hcpd, aes(x = age, y = SA.axis, group = tract, alpha = significant.derivative, color = SA.axis, linewidth = significant.derivative)) +
  scale_alpha_continuous(range = c(0, 1), breaks = c(0, 0.001, 0.002, 0.003, 0.004, 0.005, 0.006), limits = c(0.0005, 0.006)) +
  scale_linewidth_continuous(range = c(0, 4.5), breaks = c(0, 0.001, 0.002, 0.003, 0.004, 0.005, 0.006),  limits = c(0.0005, 0.006)) +
  scale_color_gradient2(low = "#ffc933", mid = "#e9dcf2", high = "#6f1282", guide = "colorbar", midpoint = 180)  +
  scale_x_continuous(breaks = c(8, 10, 12, 14, 16, 18, 20, 22), expand = c(0, .25)) +
  scale_y_continuous(expand = c(0.03, 0)) +
  theme_classic() +
  ylab("Sensorimotor-Association Axis Rank\n") +
  xlab("\nAge") +
  theme(
   legend.position = "none",
   axis.line = element_line(size = .2),
   axis.text = element_text(size = 5, family = "Arial", color = c("black")),
   axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
   axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
   axis.ticks = element_line(linewidth = .2))

Thalamocortical connections mature at progressively older ages across the S-A axis

PNC

Spearman’s rho

cor.test(gameffects_FA_glasser_pnc$SA.axis, gameffects_FA_glasser_pnc$smooth.increase.offset, method = c("spearman"), exact = F)
## 
##  Spearman's rank correlation rho
## 
## data:  gameffects_FA_glasser_pnc$SA.axis and gameffects_FA_glasser_pnc$smooth.increase.offset
## S = 685321, p-value = 2.986e-13
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.4859962

Spatial permutation (spin) based p-value

spin.data <- left_join(spin.parcels, gameffects_FA_glasser_pnc, by = c("label", "orig_parcelname", "tract")) 
perm.sphere.p(x = spin.data$SA.axis, y = spin.data$smooth.increase.offset, perm.id = perm.id.full, corr.type = "spearman") 
## [1] 8e-04

Correlation plot

ggplot(gameffects_FA_glasser_pnc, aes(x = SA.axis, y = smooth.increase.offset, color = SA.axis)) +
  geom_point(size = 0.5) +
  geom_smooth(method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) + 
  scale_color_gradient2(low = "goldenrod1", mid = "#ede4f5", high = "#6f1282", guide = "colorbar", aesthetics = "color", midpoint = 180) +
  theme_classic() +
  scale_y_continuous(limits = c(12, 23), breaks = c(12, 14, 16, 18, 20, 22)) +
  labs(x="\nS-A axis", y="Age of maturation (years)\n") +
  theme(
  legend.position = "none", 
  axis.text = element_text(size = 5, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2)) 

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure5/SAaxis_agematuration_pnc.pdf", device = "pdf", dpi = 500, width = 1.93, height = 1.62)

Brain map

Thalamocortical Age of Maturation

agemat.df <- gameffects_FA_glasser_pnc %>% select(label, smooth.increase.offset)
agemat.df <- left_join((spin.data %>% select(label)), agemat.df, by = "label")
agemat.df$cortex <- "cortex"
agemat.df <- rbind(agemat.df, data.frame(label = "rh_???", smooth.increase.offset = NA, cortex = "medialwall")) #create a label for medial wall to color it grey
agemat.df <- rbind(agemat.df, data.frame(label = "lh_???", smooth.increase.offset = NA, cortex = "medialwall")) #create a label for medial wall to color it grey

ggplot() + 
  geom_brain(data = agemat.df, atlas = glasser, mapping = aes(fill = smooth.increase.offset, colour=I("black"), size=I(.05))) +  
  theme_void()  + 
  scale_color_gradient2(low = "goldenrod1", mid = "#c4b0d6", high = "#6f1282", guide = "colorbar", aesthetics = "fill", midpoint =  17.5, na.value = "gray93", limits = c(16, 19), oob = squish) + 
  new_scale_fill() + 
  geom_brain(data = agemat.df, atlas = glasser, mapping = aes(fill = cortex, colour=I("black"), size=I(.05))) + 
  scale_fill_manual(values = c(alpha("white", 0), "white")) +
  theme(
   legend.text = element_text(size = 5, family = "Arial", color = c("black")),
   legend.key.size = unit(1, 'mm'),
   legend.title = element_text(size = 5, family = "Arial", color = c("black")))

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure5/Ageofmaturation_corticalplot_pnc.pdf", device = "pdf", dpi = 500, width = 5, height = 6.5) 

Sensorimotor-Association Axis

S.A.axis.plot <- S.A.axis[S.A.axis$tract %in% tractlist.pnc$tract,] %>% select(label, SA.axis)

ggplot() + 
  geom_brain(data = S.A.axis.plot, atlas = glasser, mapping = aes(fill = SA.axis, colour=I("black"), size=I(.05))) +  
  theme_void()  + 
  scale_color_gradient2(low = "goldenrod1", mid = "#c4b0d6", high = "#6f1282", guide = "colorbar", aesthetics = "fill", midpoint =  180, na.value = "gray93") + 
   new_scale_fill() + 
  geom_brain(data = agemat.df, atlas = glasser, mapping = aes(fill = cortex, colour=I("black"), size=I(.05))) + 
  scale_fill_manual(values = c(alpha("white", 0), "white")) +
  theme(
   legend.text = element_text(size = 5, family = "Arial", color = c("black")),
   legend.key.size = unit(1, 'mm'),
   legend.title = element_text(size = 5, family = "Arial", color = c("black")))

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure5/SAaxis_corticalplot_pnc.pdf", device = "pdf", dpi = 500, width = 4.7, height = 6.2) 

HCPD

Spearman’s rho

cor.test(gameffects_FA_glasser_hcpd$SA.axis, gameffects_FA_glasser_hcpd$smooth.increase.offset, method = c("spearman"), exact = F)
## 
##  Spearman's rank correlation rho
## 
## data:  gameffects_FA_glasser_hcpd$SA.axis and gameffects_FA_glasser_hcpd$smooth.increase.offset
## S = 342656, p-value = 2.008e-12
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.5164053

Spatial permutation (spin) based p-value

spin.data <- left_join(spin.parcels, gameffects_FA_glasser_hcpd, by = c("label", "orig_parcelname", "tract")) 
perm.sphere.p(x = spin.data$SA.axis, y = spin.data$smooth.increase.offset, perm.id = perm.id.full, corr.type = "spearman") 
## [1] 0.00185

Correlation plot

ggplot(gameffects_FA_glasser_hcpd, aes(x = SA.axis, y = smooth.increase.offset, color = SA.axis)) +
  geom_point(size = 0.5) +
  geom_smooth(method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) + 
  scale_color_gradient2(low = "goldenrod1", mid = "#ece4f2", high = "#6f1282", guide = "colorbar", aesthetics = "color", midpoint = 180) +
  theme_classic() +
  scale_y_continuous(limits = c(10.7, 22), breaks = c(12, 14, 16, 18, 20, 22)) +
  labs(x="\nS-A axis", y="Age of maturation (years)\n") +
  theme(
  legend.position = "none", 
  axis.text = element_text(size = 5, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2)) 

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure5/SAaxis_agematuration_hcpd.pdf", device = "pdf", dpi = 500, width = 1.93, height = 1.62)

S-A axis age resolved analysis

PNC

#Calculate the median derivative - S-A axis correlation value and its 95% credible interval at 100 ages 
colnames(temporalSAcorr_FA_glasser_pnc) <- sprintf("draw%s", seq(from = 1, to = ncol(temporalSAcorr_FA_glasser_pnc))) #set column names; 100 age x 10,000 posterior draw correlations df

deriv.SAaxis.mediancorr <- apply(temporalSAcorr_FA_glasser_pnc, 1, function(x){median(x)}) #median correlation value at each of the 100 ages
cor.credible.intervals <- apply(temporalSAcorr_FA_glasser_pnc, 1, function(x){quantile(x, probs = c(0.025, 0.975))}) #compute the credible interval for the correlation value at each age based on all draws

cor.credible.intervals <- t(cor.credible.intervals)
cor.credible.intervals <- as.data.frame(cor.credible.intervals)
cor.credible.intervals$age <- as.numeric(seq(8, 23, length.out = 100))
cor.credible.intervals$SA.correlation <- deriv.SAaxis.mediancorr 
colnames(cor.credible.intervals) <- c("lower","upper","age","SA.correlation")
ggplot(cor.credible.intervals, aes(x = age, y = SA.correlation, ymin = lower, ymax = upper)) + 
  geom_ribbon(fill = c(alpha("#323280", 0.4))) +
  geom_line(linewidth = .2, color = "black") +
  labs(x="\nAge", y="Developmental Alignment to the S-A Axis\n") +
  theme_classic() + 
  scale_x_continuous(breaks=c(8, 12, 16, 20), expand = c(0, .1)) +
  scale_y_continuous(breaks=c(0, 0.1, 0.2, 0.3, 0.4), limits = c(0, 0.4)) +
  theme(
   axis.text = element_text(size = 5, family = "Arial", color = c("black")),
   axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
   axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
   axis.line = element_line(linewidth = .2), 
   axis.ticks = element_line(linewidth = .2)) 

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure5/Derivative-SAcorrelation_plot_pnc.pdf", device = "pdf", dpi = 500, width = 1.35, height = 1.03)

HCPD

#Calculate the median derivative - S-A axis correlation value and its 95% credible interval at 100 ages 
colnames(temporalSAcorr_FA_glasser_hcpd) <- sprintf("draw%s", seq(from = 1, to = ncol(temporalSAcorr_FA_glasser_hcpd))) #set column names; 100 age x 10,000 posterior draw correlations df

deriv.SAaxis.mediancorr <- apply(temporalSAcorr_FA_glasser_hcpd, 1, function(x){median(x)}) #median correlation value at each of the 100 ages
cor.credible.intervals <- apply(temporalSAcorr_FA_glasser_hcpd, 1, function(x){quantile(x, probs = c(0.025, 0.975))}) #compute the credible interval for the correlation value at each age based on all draws

cor.credible.intervals <- t(cor.credible.intervals)
cor.credible.intervals <- as.data.frame(cor.credible.intervals)
cor.credible.intervals$age <- as.numeric(seq(8, 23, length.out = 100))
cor.credible.intervals$SA.correlation <- deriv.SAaxis.mediancorr 
colnames(cor.credible.intervals) <- c("lower","upper","age","SA.correlation")
ggplot(cor.credible.intervals, aes(x = age, y = SA.correlation, ymin = lower, ymax = upper)) + 
  geom_ribbon(fill = c(alpha("#323280", 0.4))) +
  geom_line(linewidth = .2, color = "black") +
  labs(x="\nAge", y="Developmental Alignment to the S-A Axis\n") +
  theme_classic() + 
  scale_x_continuous(breaks=c(8, 12, 16, 20), expand = c(0, .1)) +
  #scale_y_continuous(breaks=c(0, 0.1, 0.2, 0.3, 0.4), limits = c(0, 0.4)) +
  theme(
   axis.text = element_text(size = 5, family = "Arial", color = c("black")),
   axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
   axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
   axis.line = element_line(linewidth = .2), 
   axis.ticks = element_line(linewidth = .2)) 

Thalamocortical maturation along anatomical axes

#Glasser parcel x,y,z coordinates
load(file = "/cbica/projects/thalamocortical_development/Maps/parcellations/surface/hcp_mmp1.0.rda") #RDA file from the brainGraph package https://github.com/cwatson/brainGraph/blob/master/data/hcp_mmp1.0.rda with glasser region x, y, and z MNI coordinates. regions are in lh --> rh order
hcp_mmp1.0$label <- NA
hcp_mmp1.0$label[1:180] <- glasser.regions$label[181:360]
hcp_mmp1.0$label[181:360] <- glasser.regions$label[1:180]
#Function to compute the correlation between a development map and an anatomical axis (x,y,z) as well as test the significance of the difference between the developmental correlation with the anatomical axis versus the S-A axis
compute_axis_correlation <- function(df.dev, dev.metric, axis){
  df.dev.axes <- merge(df.dev, hcp_mmp1.0, by = "label")
  df.dev.axes.spin <- left_join(spin.parcels, df.dev.axes,  by = c("label", "tract", "orig_parcelname"))
  
  #S-A axis - development correlation (as computed above) for comparison of two overlapping correlations based on dependent groups
  SA.axis.cor <-  cor(df.dev.axes$SA.axis, df.dev.axes[,dev.metric], use = "complete.obs", method = c("spearman"))
  
  #Anatomical axis - development correlation
  anatomical.axis.cor <- cor(df.dev.axes[,axis], df.dev.axes[,dev.metric], use = "complete.obs", method = c("spearman")) #correlation between anatomical axis and development metric
  anatomical.axis.pvalue <- perm.sphere.p(x = df.dev.axes.spin[,axis], y = df.dev.axes.spin[,dev.metric], perm.id = perm.id.full, corr.type = "spearman") #spin based p-value for the correlation
  
  #S-A axis - anatomical axis correlation
  SA.anatomical.cor <- cor(df.dev.axes$SA.axis, df.dev.axes[,axis], use = "complete.obs", method = c("spearman"))
  
  #Test for significant difference between two correlations based on dependent groups (i.e., correlations with an overlapping input)
  cocor.result <- cocor.pvalue <- cocor.dep.groups.overlap(
  r.jk = SA.axis.cor, 
  r.jh = anatomical.axis.cor, 
  r.kh = SA.anatomical.cor, 
  n = nrow(df.dev.axes), alternative = "two.sided", test = "hittner2003")
  
  cocor.pvalue <- cocor.result@hittner2003[["p.value"]]
  
  comparison.results <- data.frame(axis = axis, axis.cor = anatomical.axis.cor, axis.cor.pvalue = anatomical.axis.pvalue, SA.cor = SA.axis.cor, cocor.pvalue = cocor.pvalue)
  
  return(comparison.results)
  }

PNC

Anterior-posterior axis

compute_axis_correlation(df.dev = gameffects_FA_glasser_pnc, dev.metric = "smooth.increase.offset", axis = "y.mni")
##    axis  axis.cor axis.cor.pvalue    SA.cor cocor.pvalue
## 1 y.mni 0.3097694         0.06185 0.4859962 0.0001762558

Medial-lateral axis

compute_axis_correlation(df.dev = gameffects_FA_glasser_pnc, dev.metric = "smooth.increase.offset", axis = "x.mni")
##    axis    axis.cor axis.cor.pvalue    SA.cor cocor.pvalue
## 1 x.mni -0.09059357         0.56025 0.4859962 1.604827e-10

Dorsal-ventral axis

compute_axis_correlation(df.dev = gameffects_FA_glasser_pnc, dev.metric = "smooth.increase.offset", axis = "z.mni")
##    axis  axis.cor axis.cor.pvalue    SA.cor cocor.pvalue
## 1 z.mni 0.1192805          0.2582 0.4859962  0.000236177

FDR-corrected p-values across cocor comparisons

anatomical.ps <- c(0.0001762558, 1.604827e-10, 0.000236177)
p.adjust(anatomical.ps, method = c("fdr"))
## [1] 2.361770e-04 4.814481e-10 2.361770e-04
#Put above results into df for plotting
axis.results <- data.frame(Axis = factor(), corr = double())
axis.results <- axis.results %>% add_row(Axis = "S-A", corr = 0.4859962)
axis.results <- axis.results %>% add_row(Axis = "A-P", corr = 0.3097694)
axis.results <- axis.results %>% add_row(Axis = "D-V", corr = 0.1192805)
axis.results <- axis.results %>% add_row(Axis = "M-L", corr = -0.09059357)
axis.results$Axis <- factor(axis.results$Axis, ordered = TRUE, levels = c("S-A", "A-P", "D-V", "M-L"))

ggplot(axis.results, aes(x = Axis, y = corr, fill = Axis)) + 
  geom_bar(stat='identity') + 
  labs(x="") +
  labs(y="\nAge of Maturation Correlation") +
  theme_classic()+
  scale_fill_manual(values=c(alpha("#323280", 0.5), "#672975", "#672975", "#672975")) +
  theme(
  legend.position = "none", 
  axis.text = element_text(size = 5, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2))

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure5/Anatomicalaxes_agematuration_comparison_pnc.pdf", device = "pdf", dpi = 500, width = 1.4, height = 1)

HCPD

Anterior-posterior axis

compute_axis_correlation(df.dev = gameffects_FA_glasser_hcpd, dev.metric = "smooth.increase.offset", axis = "y.mni")
##    axis axis.cor axis.cor.pvalue    SA.cor cocor.pvalue
## 1 y.mni  0.47926         0.01445 0.5164053    0.3912361

Medial-lateral axis

compute_axis_correlation(df.dev = gameffects_FA_glasser_hcpd, dev.metric = "smooth.increase.offset", axis = "x.mni")
##    axis    axis.cor axis.cor.pvalue    SA.cor cocor.pvalue
## 1 x.mni -0.04929832         0.72215 0.5164053 7.947443e-11

Dorsal-ventral axis

compute_axis_correlation(df.dev = gameffects_FA_glasser_hcpd, dev.metric = "smooth.increase.offset", axis = "z.mni")
##    axis   axis.cor axis.cor.pvalue    SA.cor cocor.pvalue
## 1 z.mni 0.06423155          0.4144 0.5164053 4.340955e-06

FDR-corrected p-values across cocor comparisons

anatomical.ps <- c(0.3912361, 7.947443e-11, 4.340955e-06)
p.adjust(anatomical.ps, method = c("fdr"))
## [1] 3.912361e-01 2.384233e-10 6.511432e-06

Thalamocortical Connectivity Variation Developmentally Increases in Association Cortex

#Function to compute the coefficient of variation across thalamocortical connections within bins of the S-A axis
compute_SAbin_cov <- function(df.dwi, tractlist){
  df.dwi.long <- df.dwi %>% melt(., id.vars=c("rbcid", "age", "sex", "mean_fd"),  variable.name = "tract", value.name = "dwi.metric") %>% merge(., S.A.axis, by = "tract", sort = F) 
  df.dwi.long <- df.dwi.long[df.dwi.long$tract %in% tractlist$tract,] #only consider tracts in this dataset's final tractlist
  
  cov.SAbin <- df.dwi.long %>% group_by(rbcid, SA.axis.bin) %>% #for every individual and every bin of the S-A axis
    do(cov = (cv(as.numeric(.$dwi.metric), na.rm = T))) %>% #calculate the coefficient of variation for the dwi metric within each bin
    unnest(cols = c(cov))
  cov.SAbin <- merge(cov.SAbin, (df.dwi %>% select(rbcid, age, sex, mean_fd)), by = "rbcid")
  
  return(cov.SAbin) #df with rows length of dataset N*5 and rows of rbcid, SA.axis.bin, cov
}

PNC

cov.SAbins.pnc <- compute_SAbin_cov(df.dwi = FA.glasser.pnc, tractlist = tractlist.pnc)
#gam interaction test to determine whether changes in cov by age vary across S-A axis bins
summary(gam(cov ~ s(age, k = 4) + s(SA.axis.bin, k = 4) + ti(age, SA.axis.bin, k = 4), method = "REML", data = cov.SAbins.pnc))
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## cov ~ s(age, k = 4) + s(SA.axis.bin, k = 4) + ti(age, SA.axis.bin, 
##     k = 4)
## 
## Parametric coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 0.1315122  0.0002498   526.4   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                       edf Ref.df       F  p-value    
## s(age)              1.001  1.002   12.68 0.000371 ***
## s(SA.axis.bin)      2.998  3.000 3191.64  < 2e-16 ***
## ti(age,SA.axis.bin) 3.116  3.585   13.74  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.627   Deviance explained = 62.8%
## -REML = -14561  Scale est. = 0.00035727  n = 5725
gamresults.bins <- vector("list", length = 5)
for(bin in 1:5){
  cov.SAbin.pnc <- cov.SAbins.pnc %>% filter(SA.axis.bin == bin)
  gamresults.bins[[bin]] <- gam.fit.smooth(measure = "cov", atlas = "SAbin", dataset = "pnc", region = "cov", smooth_var = "age", covariates = "sex + mean_fd", knots = 3, set_fx = TRUE) %>% select(-tract) %>% mutate(SA.bin = bin)
}
gamresults.cov.pnc <- do.call(rbind, gamresults.bins)

gamresults.cov.pnc <- gamresults.cov.pnc %>% mutate(fdr.pvalue = p.adjust(Anova.smooth.pvalue, method = c("fdr")))

knitr::kable(gamresults.cov.pnc)
GAM.smooth.Fvalue GAM.smooth.pvalue GAM.smooth.partialR2 Anova.smooth.pvalue smooth.change.onset smooth.peak.change smooth.decrease.onset smooth.decrease.offset smooth.increase.onset smooth.increase.offset smooth.last.change SA.bin fdr.pvalue
8.8422351 0.0001546 -0.0152757 0.0001546 8.166667 8.576633 8.166667 16.06784 NA NA 15.99330 1 0.0003866
0.0355421 0.9650831 0.0000624 0.9650831 NA NA NA NA NA NA NA 2 0.9650831
8.4049572 0.0002379 0.0145313 0.0002379 11.371859 11.483668 NA NA 11.371859 16.51508 16.44054 3 0.0003965
60.9790643 0.0000000 0.0966420 0.0000000 8.166667 8.762982 NA NA 8.166667 19.19849 19.12395 4 0.0000000
3.5096789 0.0302301 0.0061197 0.0302301 14.651591 15.844221 NA NA 14.651591 15.91876 15.84422 5 0.0377876
#Function to plot participant-level cov by age in a specific bin of the S-A axis
plot_SAbin_cov <- function(SA.bin, plotcolor){
  cov.SAbin.pnc <- cov.SAbins.pnc %>% filter(SA.axis.bin == SA.bin)
  
 cov_plot <- ggplot(cov.SAbin.pnc, aes(x = age, y = cov)) + 
   geom_point(color = alpha(plotcolor, 0.4), size = 0.5) +
   geom_smooth(method = "gam", formula = y ~ s(x, k = 3, fx = T), linewidth = 0.2, color = "black") +
   theme_classic() +
   labs(x="\nAge", y="Across-tract variation in FA (cov)\n") +
   scale_x_continuous(breaks = c(8, 12, 16, 20), expand = c(0,.45)) +
   scale_y_continuous(breaks = c(0.04, 0.08, 0.12, 0.16, 0.20, 0.24)) +
   theme(
    legend.position = "none", 
    axis.text = element_text(size = 5, family = "Arial", color = c("black")),
    axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
    axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
    axis.line = element_line(linewidth = .2), 
    axis.ticks = element_line(linewidth = .2)) 

  return(cov_plot)
 }
plot_SAbin_cov(SA.bin = 1, plotcolor = "goldenrod1")

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure5/SAbin1_cov_development_pnc.pdf", device = "pdf", dpi = 500, width = 1.67, height = 1.3)
plot_SAbin_cov(SA.bin = 2, plotcolor = "#FFD16A")

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure5/SAbin2_cov_development_pnc.pdf", device = "pdf", dpi = 500, width = 1.67, height = 1.3)
plot_SAbin_cov(SA.bin = 3, plotcolor = "#D7BCDA")

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure5/SAbin3_cov_development_pnc.pdf", device = "pdf", dpi = 500, width = 1.67, height = 1.3)
plot_SAbin_cov(SA.bin = 4, plotcolor ="#9859A4")

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure5/SAbin4_cov_development_pnc.pdf", device = "pdf", dpi = 500, width = 1.67, height = 1.3)
plot_SAbin_cov(SA.bin = 5, plotcolor = "#6f1282")

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure5/SAbin5_cov_development_pnc.pdf", device = "pdf", dpi = 500, width = 1.67, height = 1.3)

HCPD

cov.SAbins.hcpd <- compute_SAbin_cov(df.dwi = FA.glasser.hcpd, tractlist = tractlist.hcpd)
#gam interaction test to determine whether changes in cov by age vary across S-A axis bins
summary(gam(cov ~ s(age, k = 4) + s(SA.axis.bin, k = 4) + ti(age, SA.axis.bin), method = "REML", data = cov.SAbins.hcpd))
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## cov ~ s(age, k = 4) + s(SA.axis.bin, k = 4) + ti(age, SA.axis.bin)
## 
## Parametric coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 0.1139555  0.0002661   428.3   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                       edf Ref.df        F p-value    
## s(age)              1.404  1.685    6.056 0.00296 ** 
## s(SA.axis.bin)      2.996  3.000 1895.366 < 2e-16 ***
## ti(age,SA.axis.bin) 6.540  8.427    6.289 < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.668   Deviance explained = 66.9%
## -REML = -8065.3  Scale est. = 0.0002025  n = 2860
gamresults.bins <- vector("list", length = 5)
for(bin in 1:5){
  cov.SAbin.hcpd <- cov.SAbins.hcpd %>% filter(SA.axis.bin == bin)
  gamresults.bins[[bin]] <- gam.fit.smooth(measure = "cov", atlas = "SAbin", dataset = "hcpd", region = "cov", smooth_var = "age", covariates = "sex + mean_fd", knots = 3, set_fx = TRUE) %>% select(-tract) %>% mutate(SA.bin = bin)
}
gamresults.cov.hcpd <- do.call(rbind, gamresults.bins)

gamresults.cov.hcpd <- gamresults.cov.hcpd %>% mutate(fdr.pvalue = p.adjust(Anova.smooth.pvalue, method = c("fdr")))

knitr::kable(gamresults.cov.hcpd)
GAM.smooth.Fvalue GAM.smooth.pvalue GAM.smooth.partialR2 Anova.smooth.pvalue smooth.change.onset smooth.peak.change smooth.decrease.onset smooth.decrease.offset smooth.increase.onset smooth.increase.offset smooth.last.change SA.bin fdr.pvalue
0.762425 0.4670114 -0.0026821 0.4670114 NA NA NA NA NA NA NA 1 0.4670114
0.858181 0.4244824 -0.0030180 0.4244824 NA NA NA NA NA NA NA 2 0.4670114
24.451553 0.0000000 0.0794006 0.0000000 8.083333 8.535176 NA NA 8.083333 17.12018 17.05067 3 0.0000000
28.500997 0.0000000 0.0913491 0.0000000 8.083333 8.569933 NA NA 8.083333 17.05067 16.98116 4 0.0000000
1.097592 0.3343815 0.0038566 0.3343815 NA NA NA NA NA NA NA 5 0.4670114

Thalamocortical Connectivity Maturation Mirrors Timescales of Cortical Plasticity

Plasticity marker development data

#Fluctuation amplitude age of decrease onset data from Sydnor et al., 2023, Nat Neurosci https://www.nature.com/articles/s41593-023-01282-y 
boldamplitude.agedecrease.cifti <- read_cifti("/cbica/projects/thalamocortical_development/Maps/spatiotemp_dev_plasticity/cortical_maps/AgeofDeclineOnset_FirstNegDeriv.pscalar.nii") 

boldamplitude.dev <- data.frame(orig_parcelname = names(boldamplitude.agedecrease.cifti$Parcel), amplitude.agedecline = boldamplitude.agedecrease.cifti$data)
#T1/T2 ratio rate of developmental change from Baum et al., 2022, J Neuro https://www.jneurosci.org/content/42/29/5681
myelin.maxdev.cifti <- read_cifti("/cbica/projects/thalamocortical_development/Maps/myelin_development/hcpd_n628_median_posterior_age_of_max_slope_myelination.pscalar.nii") #age of maximal T1/T2 ratio increase
myelin.roc.cifti <- read_cifti("/cbica/projects/thalamocortical_development/Maps/myelin_development/hcpd_n628_mean_posterior_annualized_roc_myelin.pscalar.nii") #annualized rate of change

myelin.dev <- data.frame(orig_parcelname = names(myelin.maxdev.cifti$Parcel), myelin.agemaxdev = myelin.maxdev.cifti$data, myelin.annualroc = myelin.roc.cifti$data)
#Generate the E:I ratio-age slope map in Glasser HCP-MMP parcels from Zhang, Larsen et al., bioRxiv, https://www.biorxiv.org/content/10.1101/2023.06.22.546023v1.full

library(ciftiTools)

#Calculate E:I age slopes in Yan100 parcels
EI.group.ages <- read.csv("/cbica/projects/thalamocortical_development/Maps/EI_development/age_Yan.csv", header = F) %>% as.data.frame %>% set_names("age") %>% mutate(age = age/12) #average age of 29 age groups
EI.group.ratio <- t(read.csv("/cbica/projects/thalamocortical_development/Maps/EI_development/EI_Yan.csv", header = F)) %>% as.data.frame() #average E:I ratio in all of the Yan100 parcels in the 29 age groups
Yan100.EI.ageslopes <- sapply(EI.group.ratio, function(x) lm(x ~ EI.group.ages$age)$coefficients[2]) %>% as.data.frame() %>% set_names("slope") #linear model coefficients for the association between age and E:I ratio in all Yan100 parcels

#Reorder Yan100 parcels to match the released version of the atlas
Yan100.region.reordering <- read.csv("/cbica/projects/thalamocortical_development/Maps/EI_development/Yanatlas_regionmapping.csv") #mapping of regions between version 
Yan100.EI.ageslopes$mapping <- Yan100.region.reordering$parcelnumber.released 
Yan100.EI.ageslopes <- Yan100.EI.ageslopes %>% arrange(mapping)

#Map parcel-level Yan100 data to fslr 32k vertex-level data
Yan100.densecifti <- read_cifti("/cbica/projects/thalamocortical_development/Maps/EI_development/100Parcels_Yeo2011_17Networks.dscalar.nii")
Yan100.lh.vertices <- Yan100.densecifti$data$cortex_left %>% as.data.frame %>% set_names("mapping") #parcel number assignments for 32492 vertices
Yan100.lh.vertices.slope <- left_join(Yan100.lh.vertices, Yan100.EI.ageslopes, by = "mapping") %>% select("slope") #age slope for all vertices
Yan100.rh.vertices <- Yan100.densecifti$data$cortex_right %>% as.data.frame %>% set_names("mapping") #parcel number assignments for 32492 vertices
Yan100.rh.vertices.slope <- left_join(Yan100.rh.vertices, Yan100.EI.ageslopes, by = "mapping") %>% select("slope") #age slope for all vertices
lh.medialwall <- read.csv("/cbica/projects/thalamocortical_development/Maps/parcellations/surface/medialwall.mask.leftcortex.csv", header = F)
rh.medialwall <- read.csv("/cbica/projects/thalamocortical_development/Maps/parcellations/surface/medialwall.mask.rightcortex.csv", header = F)
Yan100.EI.ageslopes.densecifti <- as_cifti(cortexL = Yan100.lh.vertices.slope$slope, cortexR = Yan100.rh.vertices.slope$slope, cortexL_mwall = lh.medialwall$V1, cortexR_mwall = rh.medialwall$V1) #ciftify me captain
write_cifti(Yan100.EI.ageslopes.densecifti, "/cbica/projects/thalamocortical_development/Maps/EI_development/EI_ageslopes_Yan100.dscalar.nii")

#Parcellate the fslr 32k vertex-level E:I-age slope with the glasser HCP-MMP atlas
command = "-cifti-parcellate /cbica/projects/thalamocortical_development/Maps/EI_development/EI_ageslopes_Yan100.dscalar.nii /cbica/projects/thalamocortical_development/Maps/parcellations/surface/glasser_space-fsLR_den-32k_desc-atlas.dlabel.nii COLUMN /cbica/projects/thalamocortical_development/Maps/EI_development/EI_ageslopes_glasser360.pscalar.nii -only-numeric"
ciftiTools::run_wb_cmd(command, intern = FALSE, ignore.stdout = NULL, ignore.stderr = NULL)

detach("package:ciftiTools", unload=TRUE)
EI.ageslope.cifti <- read_cifti("/cbica/projects/thalamocortical_development/Maps/EI_development/EI_ageslopes_glasser360.pscalar.nii") 
EIratio.dev <- data.frame(orig_parcelname = names(EI.ageslope.cifti$Parcel), EI.ageslope = EI.ageslope.cifti$data)
df.list <- list(boldamplitude.dev, myelin.dev, EIratio.dev) #dfs to merge
plasticity.dev <- Reduce(function(x,y) merge(x,y, all=TRUE, sort=F), df.list) 
plasticity.dev <- left_join(spin.parcels, plasticity.dev, by = "orig_parcelname")

Fluctuation amplitude decrease onset map

ggplot() + 
  geom_brain(data = plasticity.dev, atlas = glasser, mapping = aes(fill = amplitude.agedecline, colour=I("black"), size=I(.03))) +  
  theme_void()  + 
  scale_fill_gradientn(colors = c("#fcd16d", "#EA8B57","#943b90", "#581768"), limits = c(7.9, 18), oob = squish, na.value = "white") 
## merging atlas and data by 'label'
## Warning: Some data not merged properly. Check for naming errors in data:
##   atlas type  hemi  region side  label       geometry orig_parcelname
##   <chr> <chr> <chr> <chr>  <chr> <chr> <MULTIPOLYGON> <chr>          
## 1 <NA>  <NA>  <NA>  <NA>   <NA>  lh_L…          EMPTY L_10pp_ROI     
## # ℹ 5 more variables: tract <chr>, amplitude.agedecline <dbl>,
## #   myelin.agemaxdev <dbl>, myelin.annualroc <dbl>, EI.ageslope <dbl>
## Warning: Some data not merged. Check for spelling mistakes in:
##          label orig_parcelname                     tract amplitude.agedecline
## 396 lh_L_10pp      L_10pp_ROI thalamus_L_10pp_autotrack             17.18593
##     myelin.agemaxdev myelin.annualroc EI.ageslope
## 396         16.37678      0.008886078 -0.03484702

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure6/Boldamplitude_agedecline_corticalmap.pdf", device = "pdf", dpi = 500, width = 4.75, height = 5.4) 
## merging atlas and data by 'label'
## Warning: Some data not merged properly. Check for naming errors in data:
##   atlas type  hemi  region side  label       geometry orig_parcelname
##   <chr> <chr> <chr> <chr>  <chr> <chr> <MULTIPOLYGON> <chr>          
## 1 <NA>  <NA>  <NA>  <NA>   <NA>  lh_L…          EMPTY L_10pp_ROI     
## # ℹ 5 more variables: tract <chr>, amplitude.agedecline <dbl>,
## #   myelin.agemaxdev <dbl>, myelin.annualroc <dbl>, EI.ageslope <dbl>

## Warning: Some data not merged. Check for spelling mistakes in:
##          label orig_parcelname                     tract amplitude.agedecline
## 396 lh_L_10pp      L_10pp_ROI thalamus_L_10pp_autotrack             17.18593
##     myelin.agemaxdev myelin.annualroc EI.ageslope
## 396         16.37678      0.008886078 -0.03484702

T1/T2 annualized rate of increase map

ggplot() + 
  geom_brain(data = plasticity.dev, atlas = glasser, mapping = aes(fill = myelin.annualroc, colour=I("black"), size=I(.03))) +  
  theme_void()  + 
  scale_fill_gradientn(colors = c( "#581768", "#943b90", "#EA8B57", "#fcd16d"), limits = c(0.005, 0.018), oob = squish, na.value = "white") 
## merging atlas and data by 'label'
## Warning: Some data not merged properly. Check for naming errors in data:
##   atlas type  hemi  region side  label       geometry orig_parcelname
##   <chr> <chr> <chr> <chr>  <chr> <chr> <MULTIPOLYGON> <chr>          
## 1 <NA>  <NA>  <NA>  <NA>   <NA>  lh_L…          EMPTY L_10pp_ROI     
## # ℹ 5 more variables: tract <chr>, amplitude.agedecline <dbl>,
## #   myelin.agemaxdev <dbl>, myelin.annualroc <dbl>, EI.ageslope <dbl>
## Warning: Some data not merged. Check for spelling mistakes in:
##          label orig_parcelname                     tract amplitude.agedecline
## 396 lh_L_10pp      L_10pp_ROI thalamus_L_10pp_autotrack             17.18593
##     myelin.agemaxdev myelin.annualroc EI.ageslope
## 396         16.37678      0.008886078 -0.03484702

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure6/T1T2ratio_annualizedroc_corticalmap.pdf", device = "pdf", dpi = 500, width = 4.44, height = 5.13) 
## merging atlas and data by 'label'
## Warning: Some data not merged properly. Check for naming errors in data:
##   atlas type  hemi  region side  label       geometry orig_parcelname
##   <chr> <chr> <chr> <chr>  <chr> <chr> <MULTIPOLYGON> <chr>          
## 1 <NA>  <NA>  <NA>  <NA>   <NA>  lh_L…          EMPTY L_10pp_ROI     
## # ℹ 5 more variables: tract <chr>, amplitude.agedecline <dbl>,
## #   myelin.agemaxdev <dbl>, myelin.annualroc <dbl>, EI.ageslope <dbl>

## Warning: Some data not merged. Check for spelling mistakes in:
##          label orig_parcelname                     tract amplitude.agedecline
## 396 lh_L_10pp      L_10pp_ROI thalamus_L_10pp_autotrack             17.18593
##     myelin.agemaxdev myelin.annualroc EI.ageslope
## 396         16.37678      0.008886078 -0.03484702

E:I ratio magnitude of developmental decline map

ggplot() + 
  geom_brain(data = plasticity.dev, atlas = glasser, mapping = aes(fill = EI.ageslope, colour=I("black"), size=I(.03))) +  
  theme_void()  + 
  scale_fill_gradientn(colors = c("#fcd16d", "#EA8B57","#943b90", "#581768"), limits = c(-0.04, -0.0345), oob = squish, na.value = "white") 
## merging atlas and data by 'label'
## Warning: Some data not merged properly. Check for naming errors in data:
##   atlas type  hemi  region side  label       geometry orig_parcelname
##   <chr> <chr> <chr> <chr>  <chr> <chr> <MULTIPOLYGON> <chr>          
## 1 <NA>  <NA>  <NA>  <NA>   <NA>  lh_L…          EMPTY L_10pp_ROI     
## # ℹ 5 more variables: tract <chr>, amplitude.agedecline <dbl>,
## #   myelin.agemaxdev <dbl>, myelin.annualroc <dbl>, EI.ageslope <dbl>
## Warning: Some data not merged. Check for spelling mistakes in:
##          label orig_parcelname                     tract amplitude.agedecline
## 396 lh_L_10pp      L_10pp_ROI thalamus_L_10pp_autotrack             17.18593
##     myelin.agemaxdev myelin.annualroc EI.ageslope
## 396         16.37678      0.008886078 -0.03484702

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure6/EIratio_decline_corticalmap.pdf", device = "pdf", dpi = 500, width = 4.13, height = 5.08) 
## merging atlas and data by 'label'
## Warning: Some data not merged properly. Check for naming errors in data:
##   atlas type  hemi  region side  label       geometry orig_parcelname
##   <chr> <chr> <chr> <chr>  <chr> <chr> <MULTIPOLYGON> <chr>          
## 1 <NA>  <NA>  <NA>  <NA>   <NA>  lh_L…          EMPTY L_10pp_ROI     
## # ℹ 5 more variables: tract <chr>, amplitude.agedecline <dbl>,
## #   myelin.agemaxdev <dbl>, myelin.annualroc <dbl>, EI.ageslope <dbl>

## Warning: Some data not merged. Check for spelling mistakes in:
##          label orig_parcelname                     tract amplitude.agedecline
## 396 lh_L_10pp      L_10pp_ROI thalamus_L_10pp_autotrack             17.18593
##     myelin.agemaxdev myelin.annualroc EI.ageslope
## 396         16.37678      0.008886078 -0.03484702

Thalamocortical maturational patterns align to the spatiotemporal unfolding of plasticity readouts

PNC

plasticity.dev.pnc <- left_join(plasticity.dev, gameffects_FA_glasser_pnc, by = c("label", "orig_parcelname", "tract"))

Thalamocortical connectivity maturation - E:I ratio magnitude of developmental decline

Spearman’s rho

cor.test(plasticity.dev.pnc$smooth.increase.offset, plasticity.dev.pnc$EI.ageslope, method = c("spearman"), exact = F)
## 
##  Spearman's rank correlation rho
## 
## data:  plasticity.dev.pnc$smooth.increase.offset and plasticity.dev.pnc$EI.ageslope
## S = 735977, p-value = 2.893e-11
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.4480033

Spatial permutation (spin) based p-value

EI.thalamus.pnc <- perm.sphere.p(x = plasticity.dev.pnc$smooth.increase.offset, y = plasticity.dev.pnc$EI.ageslope, perm.id = perm.id.full, corr.type = "spearman") 

Correlation plot

ggplot(plasticity.dev.pnc, aes(x = EI.ageslope, y = smooth.increase.offset)) +
  geom_point(size = 0.5, color = "#fcd16d") +
  geom_smooth(method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) + 
  scale_y_continuous(limits = c(12, 23), breaks = c(12, 14, 16, 18, 20, 22)) +
  theme_classic() +
  labs(x="E/I ratio decline (age slope)\n", y="\nThalamocortical age of maturation (years)") +
  theme(
  legend.position = "none", 
  axis.text = element_text(size = 5, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2)) 

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure6/EIdecline_thalamocortical_pnc.pdf", device = "pdf", dpi = 500, width = 1.8, height = 1.35) 

Thalamocortical connectivity maturation - myelin annualized rate of growth

Spearman’s rho

cor.test(plasticity.dev.pnc$smooth.increase.offset, plasticity.dev.pnc$myelin.annualroc, method = c("spearman"), exact = F)
## 
##  Spearman's rank correlation rho
## 
## data:  plasticity.dev.pnc$smooth.increase.offset and plasticity.dev.pnc$myelin.annualroc
## S = 1927545, p-value = 3.753e-11
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##        rho 
## -0.4456946

Spatial permutation (spin) based p-value

myelin.thalamus.pnc <- perm.sphere.p(x = plasticity.dev.pnc$smooth.increase.offset, y = plasticity.dev.pnc$myelin.annualroc, perm.id = perm.id.full, corr.type = "spearman") 

Correlation plot

ggplot(plasticity.dev.pnc, aes(x = myelin.annualroc, y = smooth.increase.offset, color = myelin.annualroc)) +
  geom_point(size = 0.5, color = "#F1A45F") +
  geom_smooth(method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) + 
  scale_x_continuous(limits = c(0.005, 0.025)) +
  scale_y_continuous(limits = c(12, 23), breaks = c(12, 14, 16, 18, 20, 22)) +
  theme_classic() +
  labs(x="T1/T2 ratio annualized growth rate\n", y="\nThalamocortical age of maturation (years)", ) +
  theme(
  legend.position = "none", 
  axis.text = element_text(size = 5, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2)) 

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure6/T1T2growth_thalamocortical_pnc.pdf", device = "pdf", dpi = 500, width = 1.8, height = 1.35) 

Thalamocortical connectivity maturation - fluctuation amplitude age of decrease onset

Spearman’s rho

cor.test(plasticity.dev.pnc$smooth.increase.offset, plasticity.dev.pnc$amplitude.agedecline, method = c("spearman"), exact = F)
## 
##  Spearman's rank correlation rho
## 
## data:  plasticity.dev.pnc$smooth.increase.offset and plasticity.dev.pnc$amplitude.agedecline
## S = 710935, p-value = 2.877e-05
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##      rho 
## 0.303949

Spatial permutation (spin) based p-value

amplitude.thalamus.pnc <- perm.sphere.p(x = plasticity.dev.pnc$smooth.increase.offset, y = plasticity.dev.pnc$amplitude.agedecline, perm.id = perm.id.full, corr.type = "spearman")

Correlation plot

ggplot(plasticity.dev.pnc, aes(x = amplitude.agedecline, y = smooth.increase.offset, color = amplitude.agedecline)) +
  geom_point(size = 0.5, color = "#7B2C80") +
  geom_smooth(method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) + 
  scale_x_continuous(limits = c(8.2, 21.2), breaks = c(8, 10, 12, 14, 16, 18, 20)) +
  scale_y_continuous(limits = c(12, 23), breaks = c(12, 14, 16, 18, 20, 22)) +
  theme_classic() +
  labs(x="Amplitude decrease onset (years)\n", y="\nThalamocortical age of maturation (years)") +
  theme(
  legend.position = "none", 
  axis.text = element_text(size = 5, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2)) 

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure6/Amplitudedecline_thalamocortical_pnc.pdf", device = "pdf", dpi = 500, width = 1.8, height = 1.35) 

FDR-corrected p-values across correlations

plasticity.ps <- c(EI.thalamus.pnc, myelin.thalamus.pnc, amplitude.thalamus.pnc)
plasticity.ps
## [1] 0.00025 0.00085 0.02850
p.adjust(plasticity.ps, method = c("fdr"))
## [1] 0.000750 0.001275 0.028500

HCPD

plasticity.dev.hcpd <- left_join(plasticity.dev, gameffects_FA_glasser_hcpd, by = c("label", "orig_parcelname", "tract"))

Thalamocortical connectivity maturation - E:I ratio magnitude of developmental decline

Spearman’s rho

cor.test(plasticity.dev.hcpd$smooth.increase.offset, plasticity.dev.hcpd$EI.ageslope, method = c("spearman"), exact = F)
## 
##  Spearman's rank correlation rho
## 
## data:  plasticity.dev.hcpd$smooth.increase.offset and plasticity.dev.hcpd$EI.ageslope
## S = 388563, p-value = 1.62e-09
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.4516166

Spatial permutation (spin) based p-value

EI.thalamus.hcpd <- perm.sphere.p(x = plasticity.dev.hcpd$smooth.increase.offset, y = plasticity.dev.hcpd$EI.ageslope, perm.id = perm.id.full, corr.type = "spearman") 

Correlation plot

ggplot(plasticity.dev.hcpd, aes(x = EI.ageslope, y = smooth.increase.offset, color = EI.ageslope)) +
  geom_point(size = 0.5, color = "#fcd16d") +
  geom_smooth(method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) + 
  scale_y_continuous(limits = c(12, 23), breaks = c(12, 14, 16, 18, 20, 22)) +
  theme_classic() +
  labs(x="E/I ratio decline (age slope)\n", y="\nThalamocortical age of maturation (years)") +
  theme(
  legend.position = "none", 
  axis.text = element_text(size = 5, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2)) 

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure6/EIdecline_thalamocortical_hcpd.pdf", device = "pdf", dpi = 500, width = 1.8, height = 1.35) 

Thalamocortical connectivity maturation - myelin annualized rate of growth

Spearman’s rho

cor.test(plasticity.dev.hcpd$smooth.increase.offset, plasticity.dev.hcpd$myelin.annualroc, method = c("spearman"), exact = F)
## 
##  Spearman's rank correlation rho
## 
## data:  plasticity.dev.hcpd$smooth.increase.offset and plasticity.dev.hcpd$myelin.annualroc
## S = 1032614, p-value = 9.472e-10
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##        rho 
## -0.4573402

Spatial permutation (spin) based p-value

myelin.thalamus.hcpd <- perm.sphere.p(x = plasticity.dev.hcpd$smooth.increase.offset, y = plasticity.dev.hcpd$myelin.annualroc, perm.id = perm.id.full, corr.type = "spearman") 

Correlation plot

ggplot(plasticity.dev.hcpd, aes(x = myelin.annualroc, y = smooth.increase.offset, color = myelin.annualroc)) +
  geom_point(size = 0.5, color = "#F1A45F") +
  geom_smooth(method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) + 
  scale_x_continuous(limits = c(0.005, 0.025)) +
  scale_y_continuous(limits = c(12, 23), breaks = c(12, 14, 16, 18, 20, 22)) +
  theme_classic() +
  labs(x="T1/T2 ratio annualized growth rate\n", y="\nThalamocortical age of maturation (years)", ) +
  theme(
  legend.position = "none", 
  axis.text = element_text(size = 5, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2)) 

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure6/T1T2growth_thalamocortical_hcpd.pdf", device = "pdf", dpi = 500, width = 1.8, height = 1.35) 

Thalamocortical connectivity maturation - fluctuation amplitude age of decrease onset

Spearman’s rho

cor.test(plasticity.dev.hcpd$smooth.increase.offset, plasticity.dev.hcpd$amplitude.agedecline, method = c("spearman"), exact = F)
## 
##  Spearman's rank correlation rho
## 
## data:  plasticity.dev.hcpd$smooth.increase.offset and plasticity.dev.hcpd$amplitude.agedecline
## S = 346503, p-value = 8.766e-09
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.4416828

Spatial permutation (spin) based p-value

amplitude.thalamus.hcpd <- perm.sphere.p(x = plasticity.dev.hcpd$smooth.increase.offset, y = plasticity.dev.hcpd$amplitude.agedecline, perm.id = perm.id.full, corr.type = "spearman") 

Correlation plot

ggplot(plasticity.dev.hcpd, aes(x = amplitude.agedecline, y = smooth.increase.offset, color = amplitude.agedecline)) +
  geom_point(size = 0.5, color = "#7B2C80") +
  geom_smooth(method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) + 
  scale_x_continuous(limits = c(8.2, 21.2), breaks = c(8, 10, 12, 14, 16, 18, 20)) +
  scale_y_continuous(limits = c(12, 23), breaks = c(12, 14, 16, 18, 20, 22)) +
  theme_classic() +
  labs(x="Amplitude decrease onset (years)\n", y="\nThalamocortical age of maturation (years)") +
  theme(
  legend.position = "none", 
  axis.text = element_text(size = 5, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2)) 

ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure6/Amplitudedecline_thalamocortical_hcpd.pdf", device = "pdf", dpi = 500, width = 1.8, height = 1.35) 

FDR-corrected p-values across correlations

plasticity.ps <- c(EI.thalamus.hcpd, myelin.thalamus.hcpd, amplitude.thalamus.hcpd)
plasticity.ps
## [1] 0.00020 0.00115 0.00590
p.adjust(plasticity.ps, method = c("fdr"))
## [1] 0.000600 0.001725 0.005900